Glacier thickness and ice volume of the Northern Andes

Tropical glacier melt provides valuable water to surrounding communities, but climate change is projected to cause the demise of many of these glaciers within the coming century. Understanding the future of tropical glaciers requires a detailed record of their thicknesses and volumes, which is currently lacking in the Northern Andes. We calculate present-day (2015–2021) ice-thicknesses for all glaciers in Colombia and Ecuador using six different methods, and combine these into multi-model ensemble mean ice thickness and volume maps. We compare our results against available field-based measurements, and show that current ice volumes in Ecuador and Colombia are 2.49 ± 0.25 km3 and 1.68 ± 0.24 km3 respectively. We detected no motion on any remaining ice in Venezuela. The overall ice volume in the region, 4.17 ± 0.35 km3, is half of the previous best estimate of 8.11 km3. These data can be used to better evaluate the status and distribution of water resources, as input for models of future glacier change, and to assess regional geohazards associated with ice-clad volcanoes.

Step 1b. ice masking. We find that RGI polygons overestimate the spatial extent of the majority of glaciers in the Northern Andes, likely due to rapid glacier recession and persistent high-altitude snow cover throughout

-Calculate glacier velocities
Multipass feature tracking
the year 1 . Any ice-volume calculation method using the RGI polygons 30,35,37 is therefore likely to overestimate the true volume. Hence, we use a new ice-masking approach, which uses percentile analysis of satellite image timeseries to differentiate permanently glaciated regions from temporary snow cover. This script generates an ice index (using Sentinel-2 visible bands 2, 3, and 4, and shortwave infra-red bands 11 and 12) and a water index (using Sentinel-2 visible bands 2 and 3, and near infra-red band 8), and classifies pixels as ice or water when the respective index is below an optically calibrated threshold in 90% of individual images. This threshold filters out temporary snow cover but retains zones of persistent ice cover. Water pixels are merged with non-ice pixels and masked out, leaving a binary ice mask. Removing water pixels is important as these may otherwise be misclassified as ice: two recent global compilations mis-identify the El Altar crater lake as ice, and include it in their glacier volume calculations 35,37 . Due to its very high tephra cover, we manually delineate the Nevado del Ruiz glacier zone by hand using 3-m-resolution Planet DoveSat images collected on 10 and 12 February, 2020. We manually evaluate each glacier mask against low snow cover Sentinel-2 images from 2020 or 2021.
Step 2. Velocity map generation. We use the feature-tracking toolbox GIV to calculate 50-m-resolution glacier-surface-velocity maps for each location, using a frequency-domain multi-pass image correlator 36 . We also calculate apparent displacements over the surrounding bedrock to correct for georeferencing errors and evaluate local noise levels. We include all image pairs with a minimum temporal separation of 3 months and a maximum temporal separation of 4 years. This allows us to compile a large dataset of image pairs at each glacier, ranging from 230 pairs at Nevado del Huila to 14,198 image pairs at Sierra Nevada de Santa Marta. Calculating a large number of individual velocity maps is advantageous, as the precision of mean velocity maps improves with the number of individual velocity maps stacked 36,48 . We do not consider glacier surges, as they have not been documented in this region and the local glacier characteristics are unfavorable to their occurrence 49 .
For each Sentinel-2 image pair, we calculate displacements using iteratively reducing multipass template matching following the standard GIV workflow. We use the standard GIV reference window sizes of 400 m, 200 m, and 100 m and a 50% window overlap, for a final velocity-map resolution of 50 m. Displacements are evaluated to sub-pixel precision in the final pass using a Gaussian sub-pixel estimator 36 . We convert each pair-wise displacement map into a velocity map by dividing it by the temporal separation between images. We filter each velocity map by removing the value for pixels which meet any of the following criteria: the velocity exceeds the maximum velocity threshold of 100 m.a −1 , the signal to noise ratio is lower than 5, or the peak ratio is less than 1.3 36 . These thresholds were manually selected based on local tests to exclude the majority of pixels with erroneous velocity estimates based on comparison with neighboring pixels and external datasets 37 . We also exclude values that differ by more than 50% from their immediate neighbours (four surrounding cells) and 200% from the mean of their larger local area (25 surrounding cells), and interpolate across these now-empty pixels using the values of the remaining (i.e. valid) ice-speed pixels. We do not filter based on flow direction because most of the region's glaciers flow radially outwards from mountain peaks. To correct for possible georeferencing errors, we subtract the median velocity over non-glacierized (stable) areas from the glacier velocity in the x and y directions for each image pair. We calculate a timeseries of velocities for each pixel, and after constructing this timeseries exclude pixels having a glacier speed in excess of 1.5 standard deviations from the mean 36 . Finally, we average all individual processed velocity maps into a mean velocity map covering the entire period. We crop mean velocity maps to the updated ice mask prior to inverting for ice thickness to exclude background noise over non-glacierized terrain from ice-volume calculations.
Step 3. ice-thickness calculation. Previous intercomparisons of glacier-thickness calculation methods have shown that, while no single approach is clearly superior to all others, the average of multiple different methodologies is generally more accurate than any single method 33,50 . We therefore use an ensemble of six different methods to calculate the thickness of all glaciers in our study area. Three of these methods use glacier surface flow speeds to invert for ice thickness 32,33,36 , two of these use a basal-shear-stress-based approach 28,31,[33][34][35] , and one method uses a mass-conservation-based approach 29,30,35 . 3a. Ice-velocity-based approaches. Glacier motion occurs through a combination of internal deformation, basal sliding, and subglacial sediment deformation. Ice-surface velocities u(H) may be written as a combination of internal deformation (u d ) and basal velocity (u b ; the sum of basal sliding and subglacial sediment deformation): Here, the ice-thickness H denotes velocities (full and from internal deformation alone) that are evaluated at the ice surface. We simplify this ice-flow equation based on the characteristics of the glaciers in this study area. Field studies have not revealed extensive subglacial sediment layers, thus subglacial sediment-deformation term should be at or near zero 23,41 . Glacial sliding requires warm-based ice and can be enhanced by water pressure 51 . In the tropical Northern Andes, seasonal temperature variations are low and the majority of glacier area is located in areas with a mean annual temperature below freezing 23 . As a consequence, basal sliding u b likely accounts for a small proportion of the total glacier surface velocity. We therefore account for basal velocity u b through a correction factor, β 34 , which corresponds to the fraction of glacier motion derived from basal sliding.
As a result, we directly relate internal deformation to glacier-surface velocity, and use this to compute a closed-form relationship between (unknown) ice thickness and observed ice velocity. Ice flows under its own weight, and the rate of internal deformation is a function of the ice thickness: 32,52,53 www.nature.com/scientificdata www.nature.com/scientificdata/ Here, τ b is the basal shear stress, A c is the Arrhenius creep constant, and n is Glen's flow exponent. Our use of the basal shear stress instead of the full driving stress for glacier motion comes from the shallow-ice approximation. Through this, we assume that local stresses induced by the ice are much greater than stresses induced by lateral coupling between columns of ice. Thin ice and steep slopes are characteristic of many Andean glaciers [38][39][40][41] . Both of these enhance the dominance of the basal shear stress within the full glacier driving stress, thereby supporting our use of the shallow-ice approximation. We expand basal shear stress, τ b , into measurable parameters: Here, f is a shape factor accounting for lateral drag along the glacier margins 31,32,54 , ρ i is the ice density, g is gravitational acceleration, α is the ice-surface slope angle (averaged over a length scale long enough that longitudinal coupling along the glacier flowline becomes negligible), and H is ice-thickness. We calculate the Arrhenius creep constant based on temperature: with the constants being = . ⋅ − A 2 4 10 * c 24 , Q c = 115 kJ mol −1 , R ≈ 0.0083145 (the ideal gas constant), and T * = 273 K 53 . We combine Eqs. 2, 3 and 4 and rearrange them to solve for ice-thickness: Here, the first term contains constants and parameters and the second term contains observations obtained from GIV (u H ) and a digital elevation model (sin(α)).Thus, the only unknown required to solve for ice-thickness, H, is ice-surface velocity. We implement this equation in three different ways, one novel to this study and two based on the work of Gantayat et al. 32 . Our implementation solves Eq. 5 at each location using the full two-dimensional ice surface flow speed and topographic slope fields (VWDV model). Gantayat et al. 's approach 32 divides the glacier into 100-m elevation bands and computes mean ice-surface flow speed and topographic slope for each elevation band, which we implement in two different ways. First, we calculate elevation-band-averaged flow speed and slope for whole ice caps, defined as isolated clusters of ice-masked pixels calculated in Step 1b (G14 w model). This approach is closer to the original method 32 but may average across multiple outlet glaciers from a single ice cap. Therefore, we also use TopoToolbox 55 to divide ice caps into individual ice-drainage basins based on topographic slope. We then calculate elevation-band-averaged flow speed and slope for each elevation band within each individual ice basin (G14 b model). This approach can better honor the variable dynamics of outlet glaciers from ice caps, but may result in artificial step changes in ice thickness at the boundaries between basins.
The representative length scale over which glacier surface slope is physically significant (longitudinal coupling length) is a function of the local ice thickness. When implementing our fully distributed ice-thickness solution, we use a value of 5 times the mean ice thickness 53,56,57 . In order to solve for this without prior knowledge of ice thickness, we iterate between ice-thickness and coupling-length calculations 5 times. In tests that we ran, three iterations were always sufficient for convergence, and applying five iterations permits a factor of safety between these tests and the present application. For the Gantayat et al. 32 basin-divided and whole glacier approaches 32 , the coupling-length is accounted for by the elevation-band averaging.
3b. Basal-shear-stress-based approaches. We may rewrite Eq. 3 to produce an alternative expression relating ice-thickness directly to basal shear stress: This equation forms the basis for the 'Glacier Bed Topography' (GlabTop) approach to calculating glacier thickness 31,34,35 . An empirical relationship between glacier-surface elevation range Δz i and basal shear stress τ b is then used to compute ice thickness: with τ b measured in kPa and Δz i in km 28,34 . Similarly to the Gantayat et al. 32 basin-divided and whole glacier approaches 32 , we implement this methods on both whole ice caps (GT2 w model) and individual ice-drainage basins (GT2 b model), and iterate between ice thickness and coupling length to calculate the final ice-thickness value. (2022) 9:342 | https://doi.org/10.1038/s41597-022-01446-8 www.nature.com/scientificdata www.nature.com/scientificdata/ 3c. Mass-conservation-based approach. The conservation-of-mass approach 29 , which has been applied globally 30,35 , may be written as: with [∂h/∂t] being the change in glacier surface elevation through time, ⋅ b being the glacier surface mass balance rate, and ∇ ⋅ → q being the ice-flux divergence in the horizontal plane 29 . Integrating Eq. 8 over the entire glacier domain Ω gives: with the glacier-wide ice-flux divergence being zero. We use the apparent mass balance b 29,30 , equal to the glacier surface mass balance rate ⋅ b minus the glacier surface elevation change rate ∂h/∂t such that: ∫ Ω = .
Ω The apparent mass balance of each grid cell b i is calculated according to: acc being the vertical mass-balance gradients for the ablation area and the accumulation areas respectively 29,30,35 . z i represents the elevation of the cell and z 0 represents the elevation of the apparent equilibrium line altitude (ELA), calculated by solving Eq. 10 for z 0 (see Eq. 11) such that the glacier-integrated apparent mass balance is zero. The glacier-width-normalized mean specific ice flux, q i , is then calculated by integrating all upstream apparent mass balance measurements ( b) and dividing by the local glacier width 30 . The ice thickness, H, is then calculated as: We use the same iteration between ice-thickness and coupling-length as described in the previous methods. In order to ensure a finite and physically meaningful glacier width measurement, we do not use this method on whole ice caps, and only apply it to individual ice-drainage basins (H & F model).
Step 3d. Assessment of uncertainties. We evaluate the uncertainty in ice thickness using the Monte Carlo method, considering parameter uncertainties in Eqs. 5, 6 and 12. We conduct N = 1000 runs for each method, randomly sampling from the probability distribution of each input parameter (Table 1). Terms in Eq. 6 and the left-hand parentheses on the the right-hand side of Eqs. 5 and 12 describe ice rheology and physical parameters. To define A c , we vary temperature uniformly between 268 and 272 K based on local temperature data 12,23,38,58,59 . We keep Glen's flow exponent, n, constant at 3 53,60,61 . We vary the shape factor uniformly between 0.8 and 1: this represents a compromise between low lateral drag near the ice cap summits (while acknowledging the presence of nunataks) and higher later drag in the outlet valley glaciers 28,32,53 . We vary ice density uniformly between 743 and 917 kg/m 323,40,53,62 , which is consistent with available ice-density profiles showing a mean ice density of around 830 kg/m 3 23,40 . No estimate of sliding velocity currently exists for any glacier in the study area, so we vary the basal drag correction β within a range of 0 to 0.4, corresponding to between 0 and 40% basal sliding 23,30,53 .
For velocity-based inversions, the mean velocity for each point is derived from the process described above in Step 2, and the standard deviation is calculated from standard deviation of measured velocities over non-glaciated terrain. Equation 5 is then solved at each grid cell in the velocity matrix N times, and we extract the mean and standard-deviation ice thicknesses at each grid cell from this single-model ensemble. We repeat the same process for basal-shear-stress-based approaches using Eq. 6 and for the mass-conservation-based approach using Eq. 12, and export all mean and standard deviation maps as geotiffs.
We subdivide these ice masses by the river catchments into which they drain, and convert each into its water-equivalent volume (using a mean ice density of 873.5 kg/m 3 , and an ice density range of 830-917 kg/m 3 ) to assess the spatially variable significance of these glaciers in regulating discharge 23,40,62 . We use TopoToolbox 55 to extract a drainage network and drainage basins from the DEM, cropped to a buffered ice-mask. We define the pour points for basin delineation as the boundary of the domain, and repeat the operations for non-rectangular DEMs extending 1 km, 5 km, and 20 km beyond each ice mask. We then sum the ice volume for each drainage basin, discarding basins containing no glaciers. (2022) 9:342 | https://doi.org/10.1038/s41597-022-01446-8 www.nature.com/scientificdata www.nature.com/scientificdata/ Step 3e. Ensemble mean ice-thicknesses and ice volumes. We calculate six suites of 1000 ice-thickness maps for each of the six methods described above (6000 ice-thickness maps in total): 1. A new, fully distributed two-dimensional ice velocity inversion (VWDV model).
2. An ice-cap-wide ice velocity inversion, with elevation-band averaged surface slope and flow speed (G14 w model) 32  For each method, we calculate the mean and standard deviation from all Monte Carlo runs. We then calculate a multi-model ensemble mean ice-thickness map of each glacier as the average of the mean ice-thickness maps generated using these six methods. We provide all six individual mean ice-thickness maps alongside the multi-model ensemble mean ice-thickness map. We compare the ice-thickness maps calculated with each of our methods to two prior global compilations, F19 35 and M22 37 , which we do not include in our multi-model ensemble mean ice-thickness calculation. The accuracy and precision of each map is further discussed in the technical validation section.
For each single-method ice-thickness map and for the multi-model ensemble mean ice-thickness map, we then calculate ice volumes as a simple area-weighted sum: where H jk is either the Monte-Carlo-derived mean ice thickness at each cell or the mean of the six such means, and dx and dy are the grid resolution along each axis, which in our study remains uniform, but which for the sake of generality we include within the summation. Similarly, we calculate the standard deviation of the ice volume as an area-weighted sum of the individual standard deviations of ice thickness ( σ jk ) at each grid cell:

Data records
We compile a database of ice velocity, thickness, and volume for every glacier in the Northern Andes calculated using the six methods described above 43 . All ice-thickness maps include an additional assessment of ice-thickness uncertainty.
ice velocity. Figure 2 shows mean ice-surface velocities for selected glaciers in the Northern Andes, which range from 0 to 50 metres per year. The fastest flowing ice is located on the western flank of the Cotopaxi ice cap. Apparent velocities over non-glaciated regions, representing the local noise level, are less than 2 meters per year in all cases. This low background noise level enables us to detect motion on even very small and slow moving glaciers. We do not detect flow above background noise on two glaciers, Illiniza Sur (Ecuador) and La Corona (Venezuela), thereby preventing us from calculating meaningful ice-thickness maps and indicating a strong likelihood that only permanent snowfields remain at these locations 63,64 . Due to the small size of these glaciers (<0.1 km 2 ), they are of negligible importance to overall ice-volume calculations.
ice-thickness and volume. Figures 3 and 4 show mean ice-thicknesses maps for the same selected glaciers calculated using the new, fully-distributed velocity-based inversion and the mutli-model ensemble mean respectively. Ice volumes calculated using all six individual methods, multi-model ensemble mean ice volumes, and previous best ice volume estimates 35,37 are provided in Table 2. Nevado del Ruiz has the highest maximum and average ice-thickness. The ice cap on Volcán Cayambe has the greatest ice volume for any individual ice cap, with estimates ranging from 0.52 ± 0.05 km 3 (basal-shear-stress-based basin-divided approach 31,34 ) to 0.74 ± 0.06 km 3 (conservation-of-mass-based method 29,30 ). Three out of six methods show a greater total volume in the El Cocuy region, although this is spread across multiple distinct ice bodies. Estimates of total ice volume in Colombia range from 1.16 ± 0.05 km 3 to 1.94 ± 0.09 km 3 , with the velocityand conservation-of-mass-based ice-thickness calculations clustering near the upper bound. Estimates of total ice volume in Ecuador range from 1.90 ± 0.08 km 3 to 2.81 ± 0.10 km 3 . These estimates are around 50% lower than the previous best estimate of 3.30 km 3 for Colombia and 4.81 km 3 for Ecuador 35 ) (Fig. 5). We cannot calculate regional volumes from the most recent global assessment of glacier thickness 37 , as more than half of the www.nature.com/scientificdata www.nature.com/scientificdata/ glaciers in this region have no data coverage. The river-catchment-subdivided volumes shows that while some mountains have a near-perfect radial drainage pattern (e.g. Cotopaxi), the majority of the ice is drained eastwards for all Ecuadorian glaciers (Fig. 7)).

technical Validation
We evaluate the ice thicknesses and volumes in our database through: 1. Examination of the spatial correlation between individual ice-thickness estimates derived from different methods. 2. Comparison of ice thicknesses calculated in this study to field measurements of ice thickness. 3. Evaluation of the sensitivity of ice-volume measurements to the distribution of input parameters.
Correlation between ice-volume estimates. We evaluate the two-dimensional correlation between the six different methods used to calculate ice thicknesses (Fig. 6). The two dimensional correlation coefficient C 2D is calculated as follows: with A and B being two ice-thickness maps, and m and n being the x and y direction indices of each ice-thickness measurement. No two ice-thickness maps calculated from different methods have a correlation coefficient greater than 0.9, justifying their usage as distinct estimates of ice-thickness. The mass-conservation-based approach 29,30 and three velocity-based inversions 32,36 generally exhibit correlation coefficients greater than 0.5 between themselves. Ice-thickness maps calculated using either the basin-divided or whole-glacier basal-shear-stress-based approaches 31,34 exhibit the lowest correlation with other methods.
Comparison to field measurements. We compare our results to data from three field studies: (i) an ice-penetrating-radar-based volume estimate for Nevado del Ruiz 39-41 , (ii) an ice core drilled to bedrock on Volcán Chimborazo 38 , and (iii) multiple ice-penetrating radar lines on Volcán Antisana 42 . Overall, we note that only a small number of on-site field measurements of ice thickness and volume exist in the Northern Andes. In addition, comparisons between our measurements and existing data are complicated by the rapid ice loss that has www.nature.com/scientificdata www.nature.com/scientificdata/    www.nature.com/scientificdata www.nature.com/scientificdata/ occurred during the time between the field measurements and the collection of the Sentinel-2 imagery that we use in our thickness inversions.
Nevado del Ruiz glacier volumes. An ice-penetrating radar (IPR) survey of the Nevado del Ruiz ice cap was conducted in 2000 [39][40][41] . These data provide an ice-cap volume of 0.57 ± 0.20 km 3 , calculated by extrapolating  between the existing grid of IPR measurements. We calculate volumes for the Nevado del Ruiz ice cap ranging from 0.29 ± 0.07 km 3 to 0.52 ± 0.06 km 3 , with a multi-model ensemble mean ice volume of 0.45 ± 0.12 km 3 (Fig. 8). Total ice-volume loss from the Nevado del Ruiz ice cap between the years 2000 and 2015-2021, calculated from differencing a 20-year timeseries of ASTER DEMs 21 , is 0.15 ± 0.03 km 3 . Subtracting the volume loss from the 2000 IPR survey suggests a current ice volume of 0.42 ± 0.20 km 3 , within uncertainty of our remotely-sensed volume (Fig. 8a). Farinotti et al. 35  Chimborazo ice thickness. An ice core was drilled to bedrock on Chimborazo's second summit, Cumbre Veintimilla, in December 2000 38 . The measured depth to bedrock from this borehole is 54.4 m (see Fig. 8b). We calculate an ice thickness at Cumbre Veintimilla ranging from 47 ± 4 m to 66 ± 6 m, with a multi-model ensemble mean ice thickness of 53 ± 7 m. ASTER-based elevation-change measurements 21 show an ice-thickness change of +4 ± 19 m over the period 2000-2018 at Cumbre Veintimilla. This would result in a present-day ice-thickness of 58 ± 19 m at this point, matching our remotely-sensed volume (Fig. 8a). Farinotti et al. 35   www.nature.com/scientificdata www.nature.com/scientificdata/ contribute to the water supply of Ecuador's capital, Quito 5,44 (Fig. 8c,d). We compare data from our inversions at the same points as the field-based ice-thickness measurements.
For outlet glacier 12, we find a mean ice thickness across all points where IPR was collected ranging from 46 ± 8 m to 60 ± 9 m. Our velocity-based inversion provides a mean ice thickness of 54 ± 16 m, and the multi-model ensemble mean mean ice thickness is 54 ± 16 m at these points. All of these values are within error of the March 2014 mean ice thickness of 58 ± 18 m for all IPR measurements (Fig. 8c). Farinotti et al. 35 and Millan et al. 37 calculate a mean ice-thickness of 75 ± 7 m and 66 ± 38 m respectively for the same points. The two basal-shear-stress-based methods 31,34 , the Gantayat et al. 32 basin-divided approach 32 35 . We do not include data from Millan et al. 37 , as they did not calculate the ice thickness or volume for Chimborazo or Nevado del Ruiz 37 .
For outlet glacier 15 α, we find a mean ice thickness across all points where IPR was collected ranging from 29 ± 12 m to 41 ± 17 m. Our new velocity-based inversion provides a mean ice thickness of 34 ± 14 m, and the multi-model ensemble mean ice thickness is 37 ± 14 m at these points. All of these values are within error of the March 2014 mean ice thickness 42 of 30 ± 16 m for all IPR measurements (Fig. 8d). Farinotti et al. 35 and Millan et al. 37 calculate a mean ice thickness of 30 ± 3 m and 17 ± 14 m respectively for the same points 35,37 . All methods have a positive point-wise correlation coefficient with the IPR-derived ice thicknesses, with the two basal-shear-stress-based approaches having the lowest correlation coefficient (0.37-0.38). The mass-conservation-based approach 29,30 and our new ice-velocity inversion exhibit the highest correlation coefficients of 0.64 and 0.55 respectively, and the multi-model ensemble mean ice thickness has a correlation coefficient of 0. 43.
The comparison presented here shows that our remotely sensed ice thicknesses and volumes are consistent with field-based data. Certain models provide a better match to the data than others, with the multi-model ensemble mean, mass-conservation-based approach 29,30 , and distributed ice-velocity inversion showing the best results. All approaches requiring division of ice caps into individual glacier basins result in artificial step changes in thickness at the boundaries between basins. We therefore recommend that studies requiring the two-dimensional pattern of glacier thickness use the ice-thickness maps produced using our new fully-distributed velocity-based inversion, while studies requiring only the glacier-wide or regional ice volume use either our new fully-distributed velocity-based inversion or the multi-model ensemble mean ice thickness maps 33,65 . Sensitivity analysis. To further evaluate potential sources of bias within our ice-thickness dataset, we run a comprehensive sensitivity analysis on all parameters used for the ice-thickness calculation. For each parameter evaluated, we hold all other parameters constant at their mean value, vary the chosen parameter within the range used in this study (Table 1), and evaluate the sensitivity of glacier-integrated ice volume for all parameters. We vary temperature (268 to 272 K), ice density (743 to 917 kg/m 3 ), the lateral drag factor (0.8 to 1), the basal sliding  Table 1. For the purpose of this sensitivity test, we vary the value of Glen's flow exponent, though in the remainder of our investigation we hold it at a constant value of 3. Comp = multimodel ensemble mean ice-thickness; H F = mass-conservation-based approach 29,30 ; VWDV = fully distributed velocity-based inversion from this study; GT2 b = basal-shear-stress-based basin-divided approach 31,34 ; GT2 w = basal-shear-stress-based whole glacier approach 31,34 ; G14 b = Gantayat et al. 32 basin-divided approach 32 ; G14 w = Gantayat et al. 32 whole glacier approach 32 . www.nature.com/scientificdata www.nature.com/scientificdata/ correction (0 to 0.4), and the vertical mass-balance gradients for the ablation area and the accumulation areas (0.008 to 0.01 and 0.004 to 0.006 respectively). We also vary the value of Glen's flow exponent from the constant value of 3 used in this study for investigative purposes, within a range of 2.9 to 3.1. For ease of comparison between different mean glacier volumes, we evaluate the sensitivity ∂H as a percentage change: with V i being the ice volume resulting from a specific parameter combination and V i(ref ) being a reference ice volume (calculated with all parameters equal to their mean value). Ice volume is most sensitive to temperature, ice density, and the lateral drag factor (f), with sensitivities of 5 to 10% to these parameters in isolation (Fig. 9). Ice density and lateral drag factor have the largest sensitivity for the basal-shear-stress-based approaches 31,34 , as they are the only parameters used by this model. Varying Glen's flow-law exponent n by 0.1 affects final ice volumes by 20% or more, highlighting the importance of this value. Changes in vertical mass-balance gradient only affect the mass-conservation-based approach 29,30 , and have a small effect (<5%). A table providing the results of the sensitivity test for each individual glacier region is available in the data record 43 .

Usage Notes
All ice-velocity, ice-thickness, ice-thickness-uncertainty, and basin-divided water-equivalent-volume grids may be downloaded from the Zenodo repository 43 . Individual ice-thickness maps are available for each glaciated area using the six individual methods described in this study, alongside the multi-model ensemble mean ice-thickness map. All data are saved in 32-bit floating-point geotiff format. The data are freely available under the Creative Commons Attribution Licence, CC BY 4.0.

Code availability
The feature-tracking code used to derive ice velocities, GIV, is available on github and Zenodo 36,66 . All other code, including Google Earth Engine download scripts and the ice-thickness inversion code, is available on zenodo 67 .